Anisotropic Parameter Determination

ABSTRACT

A walkaway VSP survey is carried out using a receiver array. First arrivals to a plurality of receivers are picked and used to estimate a normal-moveout (NMO) velocity. Using the NMO velocity and vertical velocities estimated from the VSP data, two anisotropy parameters are estimated for each of the layers. The anisotropy parameters may then be used to process surface seismic data to give a stacked image in true depth and for the interpretation purposes. For multi-azimuthal walkaway or 3D VSP data, we determine two VTI parameters ε and δ for multi-azimuth vertical planes. Then we determine five anisotropic interval parameters that describe P-wave kinematics for orthorhombic layers. These orthorhombic parameters may then be used to process surface seismic data to give a stacked image in true depth and for the interpretation purposes.

CROSS-REFERENCES TO RELATED APPLICATIONS

This application claims priority from U.S. Provisional Patent Application Ser. No. 61/117,613 filed on Nov. 25, 2008.

BACKGROUND OF THE DISCLOSURE

1. Field of the Disclosure

This disclosure relates to a method of geophysical prospecting which improves the accuracy of seismic migration and depth velocity model building. Specifically, the disclosure uses a walkaway, multi-azimuthal walkaway or 3D VSP survey for determination of anisotropy parameters characterizing subsurface velocities that may be used for imaging of reflections.

2. Description of the Related Art

In surface seismic exploration, energy imparted into the earth by a seismic source reflects from subsurface geophysical features and is recorded by a multiplicity of receivers. This process is repeated numerous times, using source and receiver configurations which may either form a line (2-D acquisition) or cover an area (3-D acquisition). The data which results is processed to produce an image of the reflectors using a procedure known as migration.

Conventional reflection seismology utilizes surface sources and receivers to detect reflections from subsurface impedance contrasts. The obtained image often suffers in spatial accuracy, resolution and coherence due to the long and complicated travel paths between source, reflector, and receiver.

Prior art methods have used a walkaway Vertical Seismic Profile (VSP) and 3D VSP surveys to estimate formation velocities. Common to the prior art methods is an assumption that the earth is isotropic. It has been recognized for several years that the earth is anisotropic, i.e., that the velocity of compressional waves depends upon the direction of propagation.

P-wave anisotropy is manifested by a change in the compressional wave velocity with direction of propagation in earth formations due to combined effects of sedimentary layering and the intrinsic anisotropy of the rock. Shales, in particular, could exhibit more than a 20% difference in P-wave velocities parallel to bedding and P-wave velocities perpendicular to bedding. Sandstones and limestones usually show smaller differences in velocity with direction of propagation. Postma (1955) showed that a type of anisotropy called transverse isotropy could be exhibited by seismic waves propagating through a thin layering of isotropic materials.

Determination of anisotropic velocities from surface seismic data using reflected waves is difficult due to the relatively poor data quality (regular and irregular noise influence when applying velocity analysis on primary waves) and the relatively low frequencies of surface seismic data. This procedure requires long source-receiver offsets, about 1.5 or more times the depth of interest. The long offsets complicate acquisition and processing. Nevertheless, there is prior art on the determination of an anisotropic velocity model for depth imaging of seismic data. See, for example, U.S. Pat. No. 6,864,890 to Meek et al.

The present disclosure addresses the problem of determining anisotropic formation velocities using a walkaway or 3D VSP survey. In a walkaway VSP survey, measurements are made using a plurality of receivers in a borehole responsive to excitation of one or more seismic sources at a plurality of distances from the wellbore. The estimated velocities and two VTI interval parameters ε and δ may then be used for migration of the walkaway VSP data or of surface seismic data and for interpretation purpose.

SUMMARY OF THE DISCLOSURE

One embodiment of the disclosure is a method of evaluating an earth formation. The method includes deploying an array of seismic receivers in a borehole and recording seismic signals in the array of receivers responsive to an activation of a seismic source at a plurality of positions on the surface of the earth. For a walkaway geometry, the method further includes estimating, from travel times of seismic waves from the plurality of source positions to the plurality of receivers, velocities of vertically propagating seismic waves in a plurality of intervals, and two anisotropy parameters ε and δ related to a normal moveout curve of compressional waves for the plurality of intervals. The estimated vertical velocities and the estimated interval anisotropy parameters may be used to further process surface seismic data to give a seismic image in depth for the purpose of interpretation. For a multi-azimuthal walkaway or 3D VSP geometry, the method may further include estimating, from travel times of seismic waves from the plurality of source positions around the well to the plurality of receivers, velocities of vertically propagating seismic waves in a plurality of intervals, and five orthorhombic anisotropy parameters ε₁, ε₁, δ₁, δ₂ and δ₃ related to a normal moveout velocity of compressional waves for the plurality of intervals. The estimated vertical velocities and the estimated anisotropy parameters may be used to further process surface seismic data to give a seismic image in depth for the purpose of interpretation.

Another embodiment of the disclosure is a system for evaluating an earth formation. The system includes an array of seismic receivers configured to be deployed in a borehole and record seismic signals responsive to an activation of a seismic source at a plurality of positions near the surface of the earth. The system further includes a processor configured to estimate, from travel-times of seismic waves from the plurality of source positions to the plurality of receivers, velocities of vertically propagating seismic waves in a plurality of intervals, and two anisotropy parameters ε and δ related to a normal moveout curve of compressional waves for the plurality of intervals. In case of multi-azimuthal walkaway survey or 3D VSP data, the system may further include a processor configured to estimate, from travel-times of seismic waves from the plurality of source positions around the well to the plurality of receivers, velocities of vertically propagating seismic waves in a plurality of intervals, and five orthorhombic anisotropy parameters ε₁, ε₁, δ₁, δ₂ and δ₃ related to azimuthally-dependent normal moveout curves of compressional waves for the plurality of intervals.

Another embodiment of the disclosure is a computer-readable medium accessible to at least one processor. The medium includes instructions which enable the at least one processor to estimate, using travel-times of seismic downgoing direct P waves from a plurality of source positions (walkaway geometry) to a plurality of receivers in a borehole, velocities of vertically propagating seismic P waves in a plurality of intervals, and two anisotropy parameters related to a normal moveout curve of compressional waves for the plurality of intervals. In case of 3D VSP survey, the medium may include instructions which enable the at least one processor to estimate, using travel-times of seismic downgoing direct P waves from a plurality of source positions (walkaway geometry) to a plurality of receivers in a borehole, velocities of vertically propagating seismic P waves in a plurality of intervals, and five orthorhombic anisotropy parameters related to azimuthally-dependent normal moveout curves of compressional waves for the plurality of intervals.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is best understood by reference to the attached figures in which like numerals refer to like elements, and in which:

FIG. 1 illustrates the geometry of data acquisition of a walkaway VSP according to the present disclosure;

FIG. 2 illustrates the simulated geometry of a surface seismic reflection survey using a walkaway VSP survey;

FIG. 3 shows a flow chart illustrating how the method of the present disclosure may be used for processing seismic data; and

FIG. 4 shows sectors created from 3D VSP geometry to determine VTI parameters in each sector.

DETAILED DESCRIPTION OF THE DISCLOSURE

For the present disclosure, the acquisition geometry of a walkaway VSP is illustrated in FIG. 1. Shown therein is the surface of the earth 123 with a rig 121 thereon. This may be a drilling rig or it may be a mast rig which conveys a wireline into a borehole 101. The borehole 101 penetrates layers 103, 105 . . . Positioned in the borehole 101 are seismic sensors denoted by 111 a, 111 b, 111 c, 111 d . . . Each of the sensors may include a hydrophone, a single-component geophone or a multi-component geophone. Data for a single offset VSP is typically acquired using a single seismic source such as 125 a at the surface (or within a body of water at the surface). Exemplary raypath which depicts the propagation of seismic energy from the source 125 a to the detectors 111 a and 111 d are shown.

FIG. 2 shows the simulated geometry of a surface seismic survey derived from the walkaway VSP survey of FIG. 1. A vertical borehole has been assumed, the sources have been “reflected” in the borehole to simulated receiver locations on the surface 225 a, 225 b, 225 c . . . similarly, the raypaths have also been reflected to give simulated raypaths on the left side of the borehole.

Let t(x,z) be first break time (time arrival for downgoing P wave) where x is the offset and z is a receiver depth. Then we can find the reflected traveltimes t₁(x) and t₂(x) for two virtual boundaries at the depth z₁ and z₂ of two receivers:

t ₁(x)=t(−x,z ₁)+t(x,z ₁)

t ₂(x)=t(−x,z ₂)+t(x,z ₂)

If we do not have sources at the both points −x an x (at the same distance but at the opposite sides from the well), we calculate “reflected” time t₁(x) for each pint x as: t₁(x)=2t(x,z₁).

For the reflection times, we can write the Taylor series expansion:

$\begin{matrix} {{{t^{2}(x)} = {t_{0}^{2} + \frac{x^{2}}{V_{NMO}^{2}} + {c_{2}x^{4}} + \ldots}}\mspace{14mu},} & (1) \end{matrix}$

where t(x) is the time at a source receiver offset of x, (note that this distance is twice the value of x for the VSP), t₀ is the zero-offset time (when the source coincides with the receiver at the surface), V_(NMO) is the NMO velocity. Several methods are known in the art for estimating V_(NMO). The main idea of the approximations is to represent the travel time for a layered medium by a power series in x² and to select the terms in the power series so as to match the actual travel time and the first two derivatives with respect to x². This is discussed in detail in Blias (CSEG Recorder, March 2007).

Applying the Dix relation, this gives

$\begin{matrix} {V_{Dix} = {\sqrt{\frac{{t_{0,2}V_{{NMO},2}^{2}} - {t_{0,1}V_{{NMO},1}^{2}}}{t_{0,2} - t_{0,1}}}.}} & (2) \end{matrix}$

where t_(0,1), V_(NMO,1) and t_(0,2), V_(NMO,2) correspond to the reflection boundaries comprising the horizontal layer between the depths z₁ and z₂.

-   For a vertically transversely isotropic (VTI) medium,

V _(Dix) =V _(Vert)√{square root over (1+2δ)}  (3),

where V_(Vert) is the vertical velocity and δ is an anisotropy interval parameter defined by Thomsen. Thomsen defined the parameters δ and ε as

${\delta = \frac{\left. {\left( {C_{13} + C_{44}} \right)^{2} - \left( {C_{33} - C_{44}} \right)^{2}} \right)}{2\; {C_{33}\left( {C_{33} - C_{44}} \right)}}},{ɛ = \frac{C_{11} - C_{33}}{2\; C_{33}}}$

where the elastic modulii for the TI material are given by the matrix

$\begin{bmatrix} C_{11} & {C_{11} - {2\; C_{66}}} & C_{13} & 0 & 0 & 0 \\ {C_{11} - {2\; C_{66}}} & C_{11} & C_{13} & 0 & 0 & 0 \\ C_{13} & C_{13} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{bmatrix}.$

Strictly speaking, eqn (3) is correct only when V_(Dix) is estimated using short offsets. For the purposes of this disclosure, we may refer to δ as an anisotropy parameter related to the NMO velocity and ε as an anisotropy parameter related to the horizontal velocity.

V_(Vert) in eqn (3) is measurable from the zero offset survey as the vertical velocity between depths z₁ and z₂. This then gives an estimate of the value of δ in the interval between z₁ and z₂ as:

$\begin{matrix} {\delta = {\frac{V_{Dix}^{2} - V_{Vert}^{2}}{2\; V_{Vert}^{2}}.}} & (4) \end{matrix}$

The method discussed above uses the first breaks of arrivals in the borehole.

In one embodiment of the disclosure, within the receiver array depth interval, for each pair of the receivers with some minimum distance (about 50-100 m), δ is estimated. The estimated δ is smoothed to give a smoothed δ as a function of depth. If there are boundaries in the formation within the array where the subsurface property changes significantly (e.g., lithology changes), we δ may be determined between these boundaries. Parameter δ is important because knowing δ and velocity V_(Dix), calculated thorough Dix formula from surface seismic data, it is possible to estimate the vertical velocity that is needed for time-to-depth conversion of seismic post-stack data. This is discussed with reference to FIG. 3.

Walkaway VSP data are recorded 301 as discussed above. From the walkaway VSP data, V_(NMO) is estimated for a layered model 303 using methods discussed in paragraph [0016]. From the V_(NMO) estimates, V_(DIX) is calculated 305 using eqn (2). A layered model of δ is calculated 307 using eqn. (3). This may be repeated for additional arrays in the borehole. Surface seismic data are recorded 311 at or near the location of the walkaway VSP acquisition. V_(NMO) is estimated for a plurality of depths 313. V_(DIX) is calculated 315 and using the model of δ estimated from walkaway VSP data, an estimate of V_(VERT) is derived from the surface seismic data. Those versed in the art and having benefit of the present disclosure would recognize that V_(NMO) cannot be determined for a plurality of depths from surface seismic, but only for a plurality of times. To convert times to depths, a vertical velocity or a tie of surface seismic events to VSP reflections is needed. This may not be possible. This model of V_(VERT) may then be used for time to depth conversion of stacked seismic data away from the well where the VSP data are acquired.

To estimate the Thomsen ε parameter within the receiver array interval, in one embodiment of the disclosure, the moveout function t(x) is approximated by a shifted parabola.

$\begin{matrix} {{{{t_{j}(x)}=={{t_{0,j}\left( {1 - \frac{1}{A_{j}}} \right)} + {\frac{1}{A_{j}}\sqrt{t_{0.j}^{2} + {A_{j}\frac{x^{2}}{V_{{NMO}.j}^{2}}}}\mspace{14mu} j}}} = {n - 1}},n} & (5) \end{matrix}$

using the least squares method. Here

$\begin{matrix} {{t_{0,n} = {2{\sum\limits_{k = 1}^{n}\frac{h_{k}}{v_{k}}}}},} & (6) \\ {\frac{1}{V_{{NMO},n}^{2}} = \frac{\sum\limits_{k = 1}^{n}\frac{h_{k}}{v_{P,k}}}{\sum\limits_{k = 1}^{n}{h_{k}{v_{k}\left( {1 + {2\; \delta_{k}}} \right)}}}} & (7) \\ {{A_{n} = \frac{\left( {\sum\limits_{k = 1}^{n}\frac{h_{k}}{v_{k}}} \right)\left( {\sum\limits_{k = 1}^{n}{h_{k}v_{k}^{3}B_{k}}} \right)}{\left( {\sum\limits_{k = 1}^{n}{h_{k}v_{k}}} \right)^{2}}}{{B_{k} = {1 + {8\frac{\left( {ɛ_{k} - \delta_{k}} \right)}{\left( {1 + {2\; \delta_{k}}} \right)}\left( {1 + {2\frac{\delta_{k}\gamma_{k}^{2}}{\left( {1 + {2\; \delta_{k}}} \right)\left( {1 - \gamma_{k}^{2}} \right)}}} \right)}}},}} & (8) \end{matrix}$

where n is the number of a layer between the first and second receivers, i.e., there are n layers above the second receiver, v_(k) is the vertical velocity in layer k, γ_(k) is the ratio of the shear velocity v_(s) to the compressional wave velocity in layer k.

The ratio

$\frac{\gamma^{2}}{\left( {1 - \gamma^{2}} \right)}$

lies between 0 and 0.5. For an assumed value of the

$\frac{\gamma^{2}}{\left( {1 - \gamma^{2}} \right)}$

of 0.25, this gives

$\begin{matrix} {B_{k} = {1 + {8\frac{\left( {ɛ_{k} - \delta_{k}} \right)}{\left( {1 + {2\; \delta_{k}}} \right)}\left( {1 + \frac{\delta_{k}}{2\left( {1 + {2\; \delta_{k}}} \right)}} \right)}}} & (9) \end{matrix}$

with a relative error of less than 1%.

-   For the first receiver at the depth z₁, we eqns. (6)-(9) may be     rewritten replacing n by n−1. -   After some manipulations, this gives:

$\begin{matrix} {B_{n} = {\frac{\left( {t_{0,n} - t_{0,{n - 1}}} \right)\left( {{A_{n}t_{0,n}V_{{NMO},n}^{4}} - {A_{n - 1}t_{0,{n - 1}}V_{{NMO},{n - 1}}^{4}}} \right)}{\left( {{t_{0,n}V_{{NMO},n}^{2}} - {t_{0,{n - 1}}V_{{NMO},{n - 1}}^{2}}} \right)^{2}}.}} & (10) \end{matrix}$

This then gives, using eqn (9):

$\begin{matrix} {ɛ_{n} = {\delta_{n} + {\frac{B_{n} - 1}{8\; M_{n}}.{where}}}} & (11) \\ {M_{n} = {\frac{1}{\left( {1 + {2\; \delta_{n}}} \right)}\left( {1 + \frac{\delta_{n}}{2\left( {1 + {2\; \delta_{n}}} \right)}} \right)}} & (12) \end{matrix}$

Usually, walkaway survey includes one position of a multichannel receiver array. However, if it is desired to determine interval anisotropic parameters ε and δ for the entire subsurface interval above total depth (TD), the multi-channel receiver array may be deployed at additional depths. It is then possible to determine the two anisotropic parameters ε and δ for the anisotropic depth migration of the surface seismic data for the entire subsurface above total well depth.

One embodiment of the disclosure extends the method described above to an orthorhombic horizontally layered medium. Such an orthorhombic horizontally-layered medium results from introducing, for example, vertical fractures in a VTI medium. We will use Tsvankin's notations (Tsvankin, 1997) to describe orthorhombic anisotropy. For P-wave, Tsvankin's notation includes vertical velocity and five dimensionless parameters: the VTI Thomsen's parameters δ₂ and ε₂ for [x₁, x₃] symmetry plane, δ₁ and ε₁ for [x₂, x₃] symmetry plane and parameter δ₃ for in the horizontal symmetry plane [x₁,x₃]. Knowing vertical velocity from zero-offset VSP, we can estimate anisotropic coefficients δ₁ and δ₂ and symmetry azimuth θ from NMO ellipse obtained from multi-azimuth first breaks (Blias, CSEG Recorder, May 2007).

Coefficients ε₁, ε₂, and δ₃ can be estimated through first breaks non-hyperbolic approximations assuming weak anisotropy with respect to these parameters. In this case, each vertical plane with azimuth α can be considered as VTI plane (Tsvankin, 1997). Coefficients ε(α) and δ(α) in the vertical plane with azimuth α are estimated using Dix type inversion described for VTI layered model. To estimate parameter ε₁, ε₂ and δ₃, we use equation derived by Tsvankin (1997)

ε(α_(k))=ε₁ sin⁴(α_(k)−θ)+ε₂ cos⁴(α_(k)−θ)+(2ε₂+δ₃)sin²(α_(k)−θ)cos²(α_(k)−θ)

Using the least-square method, we estimate these parameters by minimizing quadratic function ƒ(ε₁,ε₂,δ₃):

${F\left( {\delta_{1},\delta_{2},\theta_{1}} \right)} = {\sum\limits_{k = 1}^{n}\left\lbrack {{\delta \left( \alpha_{k} \right)} - {\delta_{1}{\sin^{2}\left( {\alpha_{k} - \theta} \right)}} - {\delta_{2}{\cos^{2}\left( {\alpha_{k} - \theta} \right)}}} \right\rbrack^{2}}$

To estimate ε(α) and δ(α) for several azimuths α, multi-azimuthal walkaway survey or 3D VSP data are needed. Multi-azimuthal walkaway geometry may include at least 6 walkaway lines through the well, at 30° increment. For stable estimates, it's better to have 9 walkaway lines with 20° increment. 3D VSP data provides azimuthal estimates with many sectors. Modeling shows that 12 sectors with 15 degrees each provide stable estimates of five anisotropic orthorhombic parameters ε₁, ε₂, δ₁, δ₂ and δ₃. FIG. 4 demonstrates the division the circle into 12 azimuths and added times with symmetrical offsets S₁ and S₂ to eliminate influence of velocity gradient. We also stack opposite (symmetrical) offsets to simulate surface seismic NMO curve T(x) for further Dix-type inversion for interval anisotropic parameters.

Acquisition of the data may be done using Baker Hughes' multi level receiver (MLR) that can be configured from 1 to 100 levels. For 3D VSP data, we can use 300-channel SERCEL equipment with three-component hundred level borehole tool (one vertical and two horizontal components). This greatly speeds up the data acquisition. The downhole receivers can be run in combination with other logging services, either wireline or pipe-conveyed, reducing the number of trips into the well and saving rig time. In high-angle wells, the downhole receiver can be conveyed on drill pipe or coiled tubing and also run in combination with a variety of openhole logging services greatly reducing rig time.

The present disclosure addresses the problem of determining orthorhombic anisotropic formation velocities and five orthorhombic anisotropy parameters using a 3D VSP survey. In a multi-azimuthal walkaway or 3D VSP survey, measurements are made using a plurality of receivers in a borehole responsive to excitation of one or more seismic sources at a plurality of points covering some area around the well. The estimated velocities and orthorhombic anisotropic may then be used for migration of the 3D-VSP data or of 3D surface seismic data to produce an image of the earth formation and for interpretation purpose. The produced image may be used for drilling operations in the earth formation.

The inversion and migration methodology described above may be implemented on a general purpose digital computer. As would be known to those versed in the art, instructions for the computer reside on a machine readable memory device such as ROMs, EPROMs, EAROMs, Flash Memories and Optical disks. These may be part of the computer or may be linked to the computer by suitable communication channels, and may be even at a remote location. Similarly, seismic data of the type discussed above may be stored on the computer or may be linked through suitable communication channels to the computer. The communication channels may include the Internet, enabling a user to access data from one remote location and get the instructions from another remote location to process the data. The instructions on the machine readable memory device enable the computer to access the VSP data and process the data according to the method described above.

While the foregoing disclosure is directed to the preferred embodiments of the disclosure, various modifications will be apparent to those skilled in the art. It is intended that all such variations within the scope and spirit of the appended claims be embraced by the foregoing disclosure. 

1. A method of evaluating an earth formation, the method comprising: deploying a plurality of seismic receivers in a borehole and receiving seismic signals in the array of receivers responsive to an activation of a least one seismic source at a plurality of positions on the surface of the earth; estimating, from travel times of the received seismic signals, a velocity model including velocities of vertically propagating seismic waves in a plurality of intervals, and at least two anisotropy parameters related to a normal moveout curve of compressional waves for the plurality of intervals; and using the estimated velocity model for producing an image of the earth formation.
 2. The method of claim 1 further comprising using the produced image for drilling operations.
 3. The method of claim 1 wherein the plurality of positions on the surface of the earth define a multi-azimuthal walkaway geometry, and the at least two anisotropy parameters further comprise five orthorhombic anisotropy parameters ε₁, ε₁, δ₁, δ₂ and δ₃ related to a normal moveout velocity of compressional waves for the plurality of intervals.
 4. The method of claim 1 further comprising, using, for at least one of the plurality of seismic receivers, 3-component receivers.
 5. The method of claim 1 further comprising deploying the plurality of seismic receivers on a conveyance device selected from: (i) a wireline, (ii) a pipe, and (iii) coiled tubing.
 6. The method of claim 1 further comprising estimating one of the at least two anisotropy parameters by approximating a moveout function representing the travel times by a shifted parabola.
 7. A system for evaluating an earth formation, the apparatus comprising: a plurality of seismic receivers configured to be conveyed in a borehole and receive seismic signals responsive to an activation of a least one seismic source at a plurality of positions on the surface of the earth; and at least one processor configured to: (i) estimate, from travel times of the received seismic signals, a velocity model including velocities of vertically propagating seismic waves in a plurality of intervals, and at least two anisotropy parameters related to a normal moveout curve of compressional waves for the plurality of intervals; and (ii) use the estimated velocity model for producing an image of the earth formation.
 8. The system of claim 7 wherein the at least one processor is further configured to use the produced image for drilling operations.
 9. The system of claim 7 wherein the plurality of positions on the surface of the earth define a multi-azimuthal walkaway geometry, and the at least two anisotropy parameters estimated by the at least one processor further comprise five orthorhombic anisotropy parameters ε₁, ε₁, δ₁, δ₂ and δ₃ related to a normal moveout velocity of compressional waves for the plurality of intervals.
 10. The system of claim 7 wherein at least one of the plurality of seismic receivers comprises a 3-component receiver.
 11. The system of claim 7 further comprising a conveyance device configured to convey the plurality of receivers into the borehole, the conveyance device selected from: (i) a wireline, (ii) a pipe, and (iii) coiled tubing.
 12. The system of claim 7 wherein the at least one processor is further configured to estimate the at least two anisotropy parameters by approximating a moveout function representing the travel times by a shifted parabola.
 13. A computer-readable medium product having stored thereon instructions that when read by at least one processor cause the at least one processor to execute a method, the method comprising: estimating, from travel times of seismic signals received by an array of receivers in a borehole responsive to activation of at least one seismic source at a plurality of locations on a surface of earth, a velocity model including velocities of vertically propagating seismic waves in a plurality of intervals, and at least two anisotropy parameters related to a normal moveout curve of compressional waves for the plurality of intervals; and using the estimated velocity model for producing an image of the earth formation.
 14. The computer-readable medium product of claim 13 further comprising at least one of: (i) a ROM, (ii) an EPROM, (iii) an EAROM, (iv) a flash memory, and (v) an optical disks. 